#include "AnaDef.h"
#include "BaseUtils.h"
#include "dEdxGausN.h"
#include "nevent.h"
#include "fluka.h"

#define EPSILON 1e-12

TH2D *ftof=0x0;
TH2D *fkk=0x0;
TH1D *fptSpectrum=0x0;
TH1D *f2PiEff[4] = {0x0, 0x0, 0x0, 0x0};
TH1D *fEff[4] = {0x0, 0x0, 0x0, 0x0};
TH1D *f2PiSeco[4] = {0x0, 0x0, 0x0, 0x0};
TH1D *fSeco[4] = {0x0, 0x0, 0x0, 0x0};
TF1 * fFlukaAntiproton = 0x0;
TF1 * fFlukaKminus = 0x0;

Double_t fTrigEff  = -999;
Double_t PsN = -999;
Double_t PsVrecN = -999;
Double_t PsVrecVzN = -999;

Bool_t fkSpec = kFALSE;

void ini()
{
  TString opt(dEdxGausN::Opt());
  if(!opt.Contains("MC"))
    opt+="DATA";

  if(opt.Contains("10h")){
    opt+=Form("CEN%d", dEdxGausN::GetkCEN());
  }

  printf("opt %s\n", opt.Data());

  //===
  fkSpec = kTRUE;
  if(opt.Contains("MC") || opt.Contains("FF"))
    fkSpec = kFALSE;
  //===

  printf("!!!!!!!!!!!!!!!!!!!! ini: note that normalization Nevent is only relevant for DATA, for MC or FF there should not be any spectra analyzed !!!!!!!!!!!!!!!\n");
  if(fkSpec){
  fTrigEff = 1;

  if(opt.Contains("10dMC")){
    PsN = PsNLHC10dMC;
    PsVrecN = PsVrecNLHC10dMC;
    PsVrecVzN = PsVrecVzNLHC10dMC;
    fTrigEff = 0.852;
  }
  else if(opt.Contains("10eMC")){
    PsN = PsNLHC10eMC;
    PsVrecN = PsVrecNLHC10eMC;
    PsVrecVzN = PsVrecVzNLHC10eMC;
    fTrigEff = 0.852;
  }
  else if(opt.Contains("11aMC")){
    PsN = PsNLHC11aMC;
    PsVrecN = PsVrecNLHC11aMC;
    PsVrecVzN = PsVrecVzNLHC11aMC;
    fTrigEff = 0.883;
  }
  else if(opt.Contains("10hMCCEN0")){
    PsN = PsNLHC10hMCCEN0;
    PsVrecN = PsVrecNLHC10hMCCEN0;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN0;
  }
  else if(opt.Contains("10hMCCEN1")){
    PsN = PsNLHC10hMCCEN1;
    PsVrecN = PsVrecNLHC10hMCCEN1;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN1;
  }
  else if(opt.Contains("10hMCCEN2")){
    PsN = PsNLHC10hMCCEN2;
    PsVrecN = PsVrecNLHC10hMCCEN2;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN2;
  }
  else if(opt.Contains("10hMCCEN3")){
    PsN = PsNLHC10hMCCEN3;
    PsVrecN = PsVrecNLHC10hMCCEN3;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN3;
  }
  else if(opt.Contains("10hMCCEN4")){
    PsN = PsNLHC10hMCCEN4;
    PsVrecN = PsVrecNLHC10hMCCEN4;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN4;
  }
  else if(opt.Contains("10hMCCEN5")){
    PsN = PsNLHC10hMCCEN5;
    PsVrecN = PsVrecNLHC10hMCCEN5;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN5;
  }
  else if(opt.Contains("0deDATA")){
    PsN = PsNLHC0deDATA;
    PsVrecN = PsVrecNLHC0deDATA;
    PsVrecVzN = PsVrecVzNLHC0deDATA;
    fTrigEff = 0.852;
  }
  else if(opt.Contains("10dDATA")){
    PsN = PsNLHC10dDATA;
    PsVrecN = PsVrecNLHC10dDATA;
    PsVrecVzN = PsVrecVzNLHC10dDATA;
    fTrigEff = 0.852;
  }
  else if(opt.Contains("10eDATA")){
    PsN = PsNLHC10eDATA;
    PsVrecN = PsVrecNLHC10eDATA;
    PsVrecVzN = PsVrecVzNLHC10eDATA;
    fTrigEff = 0.852;
  }
  else if(opt.Contains("11aDATA")){
    PsN = PsNLHC11aDATA;
    PsVrecN = PsVrecNLHC11aDATA;
    PsVrecVzN = PsVrecVzNLHC11aDATA;
    fTrigEff = 0.883;
  }
  else if(opt.Contains("10hDATACEN0")){
    PsN = PsNLHC10hDATACEN0;
    PsVrecN = PsVrecNLHC10hDATACEN0;
    PsVrecVzN = PsVrecVzNLHC10hDATACEN0;
  }
  else if(opt.Contains("10hDATACEN1")){
    PsN = PsNLHC10hDATACEN1;
    PsVrecN = PsVrecNLHC10hDATACEN1;
    PsVrecVzN = PsVrecVzNLHC10hDATACEN1;
  }
  else if(opt.Contains("10hDATACEN2")){
    PsN = PsNLHC10hDATACEN2;
    PsVrecN = PsVrecNLHC10hDATACEN2;
    PsVrecVzN = PsVrecVzNLHC10hDATACEN2;
  }
  else if(opt.Contains("10hDATACEN3")){
    PsN = PsNLHC10hDATACEN3;
    PsVrecN = PsVrecNLHC10hDATACEN3;
    PsVrecVzN = PsVrecVzNLHC10hDATACEN3;
  }
  else if(opt.Contains("10hDATACEN4")){
    PsN = PsNLHC10hDATACEN4;
    PsVrecN = PsVrecNLHC10hDATACEN4;
    PsVrecVzN = PsVrecVzNLHC10hDATACEN4;
  }
  else if(opt.Contains("10hDATACEN5")){
    PsN = PsNLHC10hDATACEN5;
    PsVrecN = PsVrecNLHC10hDATACEN5;
    PsVrecVzN = PsVrecVzNLHC10hDATACEN5;
  }
  else {
    printf("wrong opt %s\n", opt.Data()); exit(1);
  }
  }
}

void CorrectTOFfraction(Double_t &p0y, Double_t &e0l, Double_t &e0h,  const Double_t toty, const Double_t tofpar, const Double_t tofall)
{
  p0y = (p0y*toty+tofpar)/(toty+tofall);
  e0l = (e0l*toty)/(toty+tofall);
  e0h = (e0h*toty)/(toty+tofall); 
}

void CorrectTOFratio(Double_t &p3y, Double_t &e3l, Double_t &e3h, const Double_t p0pi,  const Double_t toty, const Double_t tofpar,const Double_t  tofpi)
{
  //only scale errors
  p3y = (p3y+tofpar/p0pi/toty)/(1+tofpi/p0pi/toty);
  e3l = e3l/(1+tofpi/p0pi/toty);
  e3h = e3h/(1+tofpi/p0pi/toty);
}

void drawpar(const Int_t ntype, const TString tn, const TString outname)
{
  TList *ll = new TList;

  ll->Add(fkk);
  if(ftof){
    ll->Add(ftof);
  }

  //correspodance between tpcip and pt
  TH1D *hnor=0, *hmpv=0, *hwid=0, *hres=0, *hchi=0;
  BaseUtils::FitSlicesY(fkk, 0x0, hnor, hmpv, hwid, hres, hchi, 10, 0, 1);//0.99, 1 are quite different!!
  ll->Add(hmpv);
  ll->Add(hwid);

  //chi2
  if(tn.Contains("fith")){
    TGraphAsymmErrors *grchi2=BaseUtils::TreeToGraph(tn,"1",(Char_t*)"xx",(Char_t*)"chidof"); ll->Add(grchi2);
  }

  //totn from horiginalyield
  TGraphAsymmErrors *gtotn=BaseUtils::TreeToGraph(tn,"1",(Char_t*)"xx",(Char_t*)"totn"); ll->Add(gtotn);

  //p0, p2, p3 =rr(ratio), p4 =cova(cov wrt pion) vs xx
  TGraphAsymmErrors *gp0[ntype];
  TGraphAsymmErrors *gp2[ntype];
  TGraphAsymmErrors *gp3[ntype];
  TGraphAsymmErrors *gp4[ntype];

  for(Int_t itype=0; itype<ntype; itype++){
    gp0[itype]=BaseUtils::TreeToGraph(tn,"1",(Char_t*)"xx",Form("p%d0", itype),   (Char_t*)"",(Char_t*)"",Form("le%d0", itype),Form("he%d0", itype),0); ll->Add(gp0[itype]);
    gp2[itype]=BaseUtils::TreeToGraph(tn,"1",(Char_t*)"xx",Form("p%d2", itype),   (Char_t*)"",(Char_t*)"",Form("le%d2", itype),Form("he%d2", itype),0); ll->Add(gp2[itype]);
    gp3[itype]=BaseUtils::TreeToGraph(tn,"1",(Char_t*)"xx",Form("rr%d", itype),   (Char_t*)"",(Char_t*)"",Form("er%d", itype), Form("er%d", itype),0); gp3[itype]->SetName(Form("gr%sp%d3", tn.Data(), itype)); ll->Add(gp3[itype]);
    gp4[itype]=BaseUtils::TreeToGraph(tn,"1",(Char_t*)"xx",Form("cova%d", itype), (Char_t*)"",(Char_t*)"",(Char_t*)"",(Char_t*)"",0); gp4[itype]->SetName(Form("gr%sp%d4", tn.Data(), itype)); ll->Add(gp4[itype]);
  }

  //p10: p0: p
  //p12: p2: p
  //p13: p3: p 
  //p14: p4: p

  //p20: p0: pt
  //p23: p3: pt =rr(ratio) * (dEtady1)/(dEtady_pion)
  //p25: FF(2)  (1/2/pi/pt) dN/dpt/dy

  //p33: p3:pt = p23/(eff-ratio-to-pi)

  //p50: fully corrected pt sepctrum

  for(Int_t itype=0; itype<ntype; itype++){
    //initiate 
    TGraphAsymmErrors * gp10 = new TGraphAsymmErrors; gp10->SetName(Form("gr%sp%d10", tn.Data(), itype));ll->Add(gp10); gp10->SetTitle("p0 : p");
    TGraphAsymmErrors * gp12 = new TGraphAsymmErrors; gp12->SetName(Form("gr%sp%d12", tn.Data(), itype));ll->Add(gp12); gp12->SetTitle("p2 : p");  

    TGraphAsymmErrors * gp13 = new TGraphAsymmErrors; gp13->SetName(Form("gr%sp%d13", tn.Data(), itype));ll->Add(gp13); gp13->SetTitle("p3 : p");
    TGraphAsymmErrors * gp14 = new TGraphAsymmErrors; gp14->SetName(Form("gr%sp%d14", tn.Data(), itype));ll->Add(gp14); gp14->SetTitle("p4 : p");  

    TGraphAsymmErrors * gp20 = new TGraphAsymmErrors; gp20->SetName(Form("gr%sp%d20", tn.Data(), itype));ll->Add(gp20); gp20->SetTitle("p0 : pt"); 
    
    TGraphAsymmErrors * gp23 = new TGraphAsymmErrors; gp23->SetName(Form("gr%sp%d23", tn.Data(), itype));ll->Add(gp23); gp23->SetTitle("p3 : pt"); 
    
    TGraphAsymmErrors * gp33 = new TGraphAsymmErrors; gp33->SetName(Form("gr%sp%d33", tn.Data(), itype));ll->Add(gp33); gp33->SetTitle("p3 : pt, corrected with efficeicny"); 

    TGraphAsymmErrors * gp50 = new TGraphAsymmErrors; gp50->SetName(Form("gr%sp%d50", tn.Data(), itype)); if(fkSpec){ll->Add(gp50);} gp50->SetTitle("fully corrected pt spectrum"); 
        

    Int_t iout=-1;
    for(Int_t rawip = 0; rawip< gp0[itype]->GetN(); rawip++){
    
      //p0
      Double_t p0x=-999, p0y=-999, e0l=-999, e0h=-999;
      gp0[itype]->GetPoint(rawip, p0x, p0y);
      e0l = gp0[itype]->GetErrorYlow(rawip);
      e0h = gp0[itype]->GetErrorYhigh(rawip);

      //skip 1/10 below xmid      
      const Int_t ix = hmpv->FindBin(p0x);
      if(dEdxGausN::FilterIX(ix)){
        continue;
      }

      //important to keep iout consecutive
      iout++;

      //p2
      Double_t p2x=-999, p2y=-999, e2l=-999, e2h=-999;
      gp2[itype]->GetPoint(rawip, p2x, p2y);
      e2l = gp2[itype]->GetErrorYlow(rawip);
      e2h = gp2[itype]->GetErrorYhigh(rawip);

      //p3
      Double_t p3x=-999, p3y=-999, e3l=-999, e3h=-999;
      gp3[itype]->GetPoint(rawip, p3x, p3y);
      e3h = gp3[itype]->GetErrorYhigh(rawip);
      e3l = e3h;

      //p4
      Double_t p4x=-999, p4y=-999, e4l=-999, e4h=-999;
      gp4[itype]->GetPoint(rawip, p4x, p4y);
      e4l = gp4[itype]->GetErrorYlow(rawip);
      e4h = gp4[itype]->GetErrorYhigh(rawip);

     
      //totn
      Double_t totx=-999, toty=-999;
      gtotn->GetPoint(rawip, totx, toty);

      //tof: correct fraction p0y, and ratio p3y
      if(ftof){
        const Double_t tofpar = ftof->GetBinContent(ix, itype+1);
        const Double_t tofpi  = ftof->GetBinContent(ix, dEdxGausN::kPION+1);
        const Double_t tofall = ftof->Integral(ix, ix, 1, 4);

        if(tofall>EPSILON){
          Double_t tmpdummy = -999, p0pi = -999;
          gp0[dEdxGausN::kPION]->GetPoint(rawip, tmpdummy, p0pi);
          if(p0pi<EPSILON){
            printf("strange p0pi<eps %e\n", p0pi);exit(1);
          }
          if( TMath::Abs(p0y/p0pi-p3y)>EPSILON ){
            printf("wrong p0pi!! %e %e %e\n", p0y, p0pi, p3y); exit(1);
          }
        
          const Double_t tmpoldp0 = p0y;
          const Double_t tmpoldp3 = p3y;
          
          //p0y
          CorrectTOFfraction(p0y, e0l, e0h,  toty, tofpar, tofall);
          
          //p3y
          CorrectTOFratio(p3y, e3l, e3h, p0pi, toty, tofpar, tofpi);
          
          printf("correcting with tof! rawip %d ix %d itype %d p0pi %e, tofpar %e tofpi %e tofall %e, toty %e, p0 %e %e, p3 %e %e\n", rawip, ix, itype, p0pi, tofpar, tofpi, tofall, toty, tmpoldp0, p0y, tmpoldp3, p3y);
        }
      }      
      
      gp10->SetPoint(iout, TMath::Power(10, p0x), p0y);
      gp10->SetPointError(iout, 0, 0, e0l, e0h);

      gp12->SetPoint(iout, TMath::Power(10, p2x), p2y);
      gp12->SetPointError(iout, 0, 0, e2l, e2h);

      gp13->SetPoint(iout, TMath::Power(10, p3x), p3y);
      gp13->SetPointError(iout, 0, 0, e3l, e3h);
     
      gp14->SetPoint(iout, TMath::Power(10, p4x), p4y);
      gp14->SetPointError(iout, 0, 0, e4l, e4h);
      
      printf("test gp10, gp12  %d %d: %d %d\n", itype, ix, gp10->GetN(), gp12->GetN());

      //=======================================================
      const Double_t tpcip = hmpv->GetBinCenter(ix);

      //check
      if(TMath::Abs(p0x-totx)>EPSILON || TMath::Abs(p0x-tpcip)>EPSILON ){
        printf("drawpar p0x != totx or p0x != tpcip %e %e %e ix %d\n", p0x, totx, tpcip, ix); exit(1);
      }
      
      if(toty<EPSILON){
        printf("drawpar toty < EPSILON!! %e\n", toty); exit(1);
      }

      //pt
      const Double_t m0 = hmpv->GetBinContent(ix);
      const Double_t w0 = hwid->GetBinContent(ix);
      const Double_t ptmean = TMath::Power(10, m0);
      //in log must use different lower and higher errors!! otherwise big bug!!
      const Double_t pte0 =  ptmean - TMath::Power(10, m0-w0);
      const Double_t pte1 =  TMath::Power(10, m0+w0) - ptmean;
      
      //p20, y: p0y
      gp20->SetPoint(iout, ptmean, p0y);
      gp20->SetPointError(iout, pte0, pte1, e0l, e0h);
      printf("test setting %s %d ( %d %d ) %e --- %e %e %e\n", gp20->GetName(), iout, rawip, ix, ptmean, p0y, e0l, e0h);

      //p23, y: p3y*etacorr
      const Double_t etacorr = dEdxGausN::RatioEtaCorrection(itype, ptmean, dEdxGausN::Eta0(), dEdxGausN::Eta1());

      const Double_t etaCorrP3y = p3y*etacorr;
      const Double_t etaCorrE3l = e3l*etacorr;
      const Double_t etaCorrE3h = e3h*etacorr;

      gp23->SetPoint(iout, ptmean, etaCorrP3y);
      gp23->SetPointError(iout, pte0, pte1, etaCorrE3l, etaCorrE3h);
      printf("test setting %s %d ( %d %d ) %e --- %e %e %e %e\n", gp23->GetName(), iout, rawip, ix, ptmean, etacorr, p3y, e3l, e3h);

      //33, y: etaCorrP3y/yEffRatio , secondary, fluka

      //note that f2PiEff is in log10 scale of pT!!
      Int_t xEffRatio = f2PiEff[itype]->FindBin(TMath::Log10(ptmean));
      while(f2PiEff[itype]->GetBinContent(xEffRatio)<EPSILON && xEffRatio>1){
        xEffRatio--;
      }
      const Double_t yEffRatio = f2PiEff[itype]->GetBinContent(xEffRatio);
      const Double_t eEffRatio = f2PiEff[itype]->GetBinError(xEffRatio);
      printf("efficiency for %d found at %d %e %e %e\n", itype,  xEffRatio, f2PiEff[itype]->GetBinCenter(xEffRatio), yEffRatio, eEffRatio);

      //seco corr = prim frac
      Int_t xSeco = fSeco[itype]->FindBin(ptmean);
      while(fSeco[itype]->GetBinContent(xSeco)<EPSILON && xSeco>1){
        xSeco--;
      }
      const Double_t seco2piY = f2PiSeco[itype]->GetBinContent(xSeco);
      const Double_t seco2piE = f2PiSeco[itype]->GetBinError(xSeco);
      if(seco2piY<EPSILON){
        printf("bad seco2piY %d %e\n", xSeco, seco2piY);exit(1);
      }

      //fluka correction for anti-proton only
      Double_t flukaCor = 1;
      if(itype == dEdxGausN::kPROTON && dEdxGausN::GetkCharge()== -1 && fFlukaAntiproton){
        flukaCor = fFlukaAntiproton->Eval(ptmean);
        if(flukaCor>1+EPSILON || flukaCor < EPSILON){
          printf("bad fluka antiproton! %e %e\n", ptmean, flukaCor);exit(1);
        }
      }
      else if(itype == dEdxGausN::kKAON && dEdxGausN::GetkCharge()== -1 && fFlukaKminus){
        flukaCor = fFlukaKminus->Eval(ptmean);
        if(flukaCor>1+EPSILON || flukaCor < EPSILON){
          printf("bad fluka kminus! %e %e\n", ptmean, flukaCor);exit(1);
        }
      }  
      printf("test setting fluka %e %e\n", ptmean, flukaCor);

      //
      const Double_t etaEffCorrP3y = etaCorrP3y/yEffRatio * seco2piY * flukaCor;
      const Double_t etaEffCorrE3l = etaEffCorrP3y*TMath::Sqrt(TMath::Power(etaCorrE3l/etaCorrP3y,2)+TMath::Power(eEffRatio/yEffRatio,2) + TMath::Power(seco2piE/seco2piY,2));
      const Double_t etaEffCorrE3h = etaEffCorrP3y*TMath::Sqrt(TMath::Power(etaCorrE3h/etaCorrP3y,2)+TMath::Power(eEffRatio/yEffRatio,2) + TMath::Power(seco2piE/seco2piY,2));

      gp33->SetPoint(iout, ptmean, etaEffCorrP3y);
      gp33->SetPointError(iout, pte0, pte1, etaEffCorrE3l, etaEffCorrE3h);
      printf("test setting %s %d ( %d %d ) %e --- %e %e %e %e -- %e %e -- %e %e\n", gp33->GetName(), iout, rawip, ix, ptmean, yEffRatio, etaEffCorrP3y, etaEffCorrE3l, etaEffCorrE3h, yEffRatio, eEffRatio, seco2piY, seco2piE);

      if(fkSpec){
      //50: ptSpectrum: normalized /= bw, i.e. dN/dpt

      const Double_t dy = 2*dEdxGausN::DeltaY(dEdxGausN::Mass(itype), ptmean, dEdxGausN::Eta0(), dEdxGausN::Eta1());
        
      const Int_t rawix = fptSpectrum->GetXaxis()->FindBin(ptmean);
        
      
      //scale by 1/Nev, 1/(2 pi pt) and  1/dy, no error propagation -> 1/N_ev (1/2 pi pt) dN/dpt/dy, and also fTrigEff * fflukaCor
      const Double_t rawfactor = 1 / (2*TMath::Pi()*ptmean) / dy / PsVrecVzN * (PsVrecN/PsN) * fTrigEff * flukaCor;
      const Double_t rawYieldY = fptSpectrum->GetBinContent(rawix) * rawfactor;
      const Double_t rawYieldE = fptSpectrum->GetBinError(rawix)   * rawfactor;

      //eff should be evaluated at PsVrecVz, namely tracking eff.
      const Double_t singleEffY = fEff[itype]->GetBinContent(xEffRatio);
      const Double_t singleEffE = fEff[itype]->GetBinError(xEffRatio);

      //secodary correction = primary fraction
      const Double_t secocorrY = fSeco[itype]->GetBinContent(xSeco);
      const Double_t secocorrE = fSeco[itype]->GetBinError(xSeco);
      if(secocorrY<1e-12){
        printf("bad seco frac %e\n", secocorrY);exit(1);
      }

      //fraction, eff
      const Double_t fullCorrP0y = rawYieldY * p0y / singleEffY * secocorrY;
      const Double_t fullCorrE0l = fullCorrP0y * TMath::Sqrt( TMath::Power(rawYieldE/rawYieldY,2) + TMath::Power(e0l/p0y,2) + TMath::Power(singleEffE/singleEffY,2) + TMath::Power(secocorrE/secocorrY,2));
      const Double_t fullCorrE0h = fullCorrP0y * TMath::Sqrt( TMath::Power(rawYieldE/rawYieldY,2) + TMath::Power(e0h/p0y,2) + TMath::Power(singleEffE/singleEffY,2) + TMath::Power(secocorrE/secocorrY,2));
      gp50->SetPoint(iout, ptmean, fullCorrP0y);
      gp50->SetPointError(iout, pte0, pte1, fullCorrE0l, fullCorrE0h);
      printf("test setting %s %d ( %d %d ) %e --- %e %e %e %e -- %e %e -- %e %e\n", gp50->GetName(), iout, rawip, ix, ptmean, rawYieldY, fullCorrP0y, fullCorrE0l, fullCorrE0h, singleEffY, singleEffE, secocorrY, secocorrE);
      }
    }
  }
  
  TFile *fout=new TFile(Form("par%s.root", outname.Data()),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();
}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("drawpar argc!=2\n"); 
    return 1;
  }
  dEdxGausN::Ini();
  dEdxGausN::HistPrint();

  ini();
  printf("fkSpec %d PsN %0.f PsVrecN %.0f PsVrecVzN %0.f fTrigEff %e\n", fkSpec, PsN, PsVrecN, PsVrecVzN, fTrigEff);
  if(fkSpec && PsN<EPSILON){
    printf("bad nevent %.0f\n", PsN);exit(1);
  }

  //=====

  //hist.root
  TFile *fileHist = new TFile("hist.root");

  fkk = (TH2D*)fileHist->Get("kk");// pt:tpcip
  if(!fkk){
    printf("drawpar error in hist.root no kk!\n"); fileHist->ls(); exit(1);
  }

  ftof = (TH2D*)fileHist->Get(Form("h%slogp", dEdxGausN::TOFTag()));
  if(ftof){
    printf("drawpar using tof!!\n");
  }

  fptSpectrum = (TH1D *)fileHist->Get("ptSpectrum");
  if(!fptSpectrum && fkSpec){
    TFile *filePtSpec = new TFile("ptSpec.root");
    if(!filePtSpec){
      printf("no ptSpectrum nor ptSpec.root!\n"); fileHist->ls(); exit(1);
    }
    fptSpectrum = (TH1D *)filePtSpec->Get("ptSpectrum");
    if(!fptSpectrum){
      printf("no ptSpectrum in ptSpec.root!\n"); filePtSpec->ls(); exit(1);
    }
  }

  const TString efffile = dEdxGausN::GetEffFile();
  TFile * fileEff = new TFile(efffile);
  if(!fileEff){
    printf("no fileEff %s\n", efffile.Data()); exit(1);
  }
  for(Int_t ii=0; ii<4; ii++){
    f2PiEff[ii] = (TH1D*) fileEff->Get(Form("Eff2Pi%d",ii));
    if(!f2PiEff[ii]){
      printf("no 2pieff %d\n", ii); fileEff->ls();exit(1);
    }

    fEff[ii] = (TH1D*) fileEff->Get(Form("eff%d",ii));
    if(!fEff[ii]){
      printf("no single eff %d\n",ii); fileEff->ls(); exit(1);
    }
  }

  const TString secofile = dEdxGausN::GetSecoFile();
  TFile * fileSeco = new TFile(secofile);
  if(!fileSeco){
    printf("no fileSeco %s\n", secofile.Data()); exit(1);
  }
  for(Int_t ii=0; ii<4; ii++){
    f2PiSeco[ii] = (TH1D*) fileSeco->Get(Form("hfracprimarypar%d2pi",ii));
    if(!f2PiSeco[ii]){
      printf("no single seco2pi %d\n",ii); fileSeco->ls(); exit(1);
    }

    fSeco[ii] = (TH1D*) fileSeco->Get(Form("hfracprimarypar%d",ii));
    if(!fSeco[ii]){
      printf("no single seco %d\n",ii); fileSeco->ls(); exit(1);
    }
  }

  if(dEdxGausN::GetkCharge()== -1 ){
    fFlukaAntiproton = TrackingEff_geantflukaCorrection(AliPID::kProton,kNegative);
    if(!fFlukaAntiproton){
      printf("no fluka proton!\n");return 1;
    }
    fFlukaKminus = TrackingEff_geantflukaCorrection(AliPID::kKaon,kNegative);
    if(!fFlukaKminus){
      printf("no fluka kminus!\n");return 1;
    }
  }

  const TString fin(argv[1]);
  TFile::Open(fin);
  TString outname(fin);
  outname.ReplaceAll("out","");
  outname.ReplaceAll("_post","");
  outname.ReplaceAll("_save","");
  outname.ReplaceAll(".root","");
  drawpar(dEdxGausN::Ntype(), fin.Contains("mn2d")?"mn2d":"fith", outname);

  return 0;
}
